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The ground-state and the transport properties of graphene subject to the potential of in-plane 
charged impurities are studied. The screening of the impurity potential is shown to be nonlinear, 
producing a fractal structure of electron and hole puddles. Statistical properties of this density 
distribution as well as the charge compressibility of the system are calculated in the leading-log 
approximation. The conductivity depends logarithmically on a, the dimensionless strength of the 
Coulomb interaction. The theory is asymptotically exact when a is small, which is the case for 
graphene on a substrate with a high dielectric constant. 



A number of recent experimental [1, 2, 3, 4, 5] and 
theoretical [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16] inves- 
tigations have studied the effect of charged impurities 
on the properties of graphene [17] at its neutrality point 
(NP). It has been shown that in response to the random 
potential of impurities electron density becomes very in- 
homogencous. A deeper physical understanding of these 
inhomogencitics is important for clarifying the true na- 
ture of disorder in this new material. Hence, an analytic 
approach to this problem is desirable. Below I present 
my results in this direction. 

As a model I assume that impurities with random 
charge ±e and total concentration rii are distributed ran- 
domly on a graphene sheet of area A. The system resides 
in a medium of dielectric constant n. The strength of the 
Coulomb interaction, U(r) = e 2 /nr, is characterized by 
the parameter a — e 2 / nhv where v is the Fermi velocity. 

My main results are as follows. For k»1 where the 
interaction is weak, the basic electronic properties 

of graphene at the NP can be computed to the leading 
order in l/£, where C 3> 1 is the solution of the equation 



£ = ln(l/(4a£)) 



(1) 



I find that the charge compressibility \o a t the NP (mea- 
sured in Ref. 5) is given by 



Xo = K/(2Tie 2 R) , R = £/{AaC) 



(2) 



where R is the screening length and £ — (2a^/rii ) _1 <C R 
has the meaning of a typical quantum uncertainty in the 
quasiparticle positions. Observables, such as the corre- 
lator S(r) = (n(r)n(O)) of the local density become non- 
trivial on scales r > £. For such r, I find 



S(r) = 
i?(r) = 



2tt£ 4 
K(r) 



M\/l-d 2 + (1 + 2i? 2 ) arcsini? 
K(£) = \ 



C. 



(3) 
(4) 



Here K{r) = (<I>(0)<j>(r)) is the correlator of the screened 
random potential I show that K(r) behaves as 



K(r) = 



HR/r) , 
2{R/r)\ 



£ < r < R, 
i?< r. 



(5) 



Equation (5) agrees with the results of Refs. 13 and 18 
in the limit a <C 1. Equation (3) is in a qualitative 
agreement with the numerical simulations [15] for a ~ 1 
(where my theory is at the border of validity) . 

Of special interest for transport and imaging experi- 
ments arc the contours n(r) = 0, which separate electron 
[n(r) > 0] and hole [n(r) < 0] regions. Only contours 
larger than some minimum size A -1 are important for 
the key observables (in imaging, A -1 is the spatial reso- 
lution; in transport, A -1 ~ £VZ, see below). To isolate 
them, we can separate the total density and the total 
potential into smoothly and rapidly varying parts, 

n(r) = n(r) + <5n(r) , $(r) = $(r) + <5$(r) , (6) 

where n(r) and $(r) contain only Fourier harmonics with 
k < A. Below I assume that l/£ < A <C l/R. I refer 
to the contours n(r) = as the p-n junctions (PNJ). I 
show that most of the PNJ are closed loops of diameter 
d ~ A -1 . They reside inside of successively larger loops, 
forming a self- similar set. The fractal dimension of 
the PNJ contours is equal to 3/2 at d < R but becomes 
7/4 at d > R. Finally, the conductivity at the NP is 

cr N p = { e 2 /h)cC, c= 0.50 ±0.05. (7) 

This result is supposed to be valid at temperatures high 
enough that weak localization effects [19] can be ne- 
glected but low enough that elastic scattering by impuri- 
ties is still the dominant current relaxation mechanism. 

A similar expression for cr^p can be deduced from 
Ref. 13. It has c = 4/-7T for a -C 1, which is larger than my 
value of c by about a factor of two. Numerical estimates 
of <tnp at a ~ 1 were reported in Ref. 16. 

Let us turn to the derivation. Let V(r) be the bare 
random potential due to impurities. In the fc-space it 
is given by Vk = U^N^, where and = 2i:ahv/k 
are the Fourier transforms of the impurity density and 

1/2 

the Coulomb potential, respectively. At < k -C n i 
we can treat them as zero-mean Gaussian random vari- 
ables with variances (N^N^r) = A/ij4,-k' and (VkVk') = 
(N k N k >). It is easy to see then that the coarse-grained 
random potential V(r) is infrared divergent [20]: 

(V 2 (r)) = ]T (VkVk,) ~ ^/2){hv/£) 2 \n{A^A) . (8) 

fe<A 
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This divergence is of course cured once the bare potential 
is replaced by a screened one, 



$k = Vk + n kf^k = Vk + (2irahv/k)n k , 
which is expected to have the following properties: 

kR 



(9) 



*k = Vk, $ k = Vk 



1 + fci? 



, 1. (10) 



The last equation follows from the definition of the ther- 
modynamic charge compressibility xo = ~ "-k/'&k and the 
relation between xo and R, Eq. (2). I wish to show that 
R is given by the second formula in Eq. (2), in particular, 
that it is a function of the impurity concentration. This 
means that the screening is nonlinear, i.e., the response 
of electrons in graphene to long-wavelength Fourier har- 
monics k < R^ 1 is non-perturbative in Vk- 

Formally, n(r) and <j>(r) can be sought by mini- 
mizing the total energy E[n] + £[n], where E[n] — 
A -1 F_knk + (l/2)C/k|^k| 2 is the electrostatic energy 
and £[n] is the sum of kinetic, exchange, and correla- 
tion energies [21]. Since the functional £[n] is nonlocal, 
nonlinear, and strictly speaking, unknown, further steps 
require approximations. We are aided by the parameter 
£: the condition £ 3> 1 ensures separation of two length 
scales: (i) the typical quantum uncertainly I over which 
the nonlocality of £ [n] is important and (ii) the screening 
length R 3> I at which the nonlinear screening sets in. As 
a result, a local density approximation £ = Jdre(n(r)) 
is valid and the total energy is minimized when 



$(r) = -e'(n(r)) . 



(11) 



The kernel s(n) can be computed by integrating out 
the rapid density fluctuations Sn(r) within the linear- 
response theory [6, 22] for a <C 1. The result is [23] 

e= 2 -^Hv\n\^ + 6e, Se = In ( (12) 



A\n\P 



A 



Here kp{r) = y/n\n(r)\ is the local Fermi momentum. 
Equation (11) can now be easily inverted to give 

i>|<l| sgn ($) |$| ._. 

n($)~ A ln ir4^ * »ft«A. (13) 

7r(fiw) 2 2P hvA 1 1 

Later, we will also need the derivative of function n(i>): 



dn 



2 m 



n(Hv) 2 2£ 2 \<f>\ 



|$| > hvA. (14) 



We have reduced the original quantum problem to the 
nonlinear integral equation (13) for the classical quantity 
n(r) [Note that <&[n] is given by Eq. (9).] I call this a 
rcnormalizcd Thomas- Fermi approximation (RTFA) [24] . 

In conventional gapped semiconductors, this would be 
about as far as one could go analytically before hav- 
ing to resort to heuristic estimates [25], computer sim- 
ulations [20], or variational methods [26]. The reason 



is as follows: the energy minimization is constrained 
by n(r) > 0, leading to a finite area fraction of deple- 
tion regions n(r) = 0. Therein Eq. (11) is replaced by 
3>(r) > [25], so that no analog of Eq. (13) exists, making 
the problem analytically intractable. In graphene, there 
is no constraint on the sign of density. Instead, the con- 
ditions for validity of Eq. (12) are fc<£ and A -C kp, 
which arc typically satisfied everywhere except near the 
PNJ (see below). However, this kind of "depiction re- 
gions" occupy a paramctrically small area fraction. 

Yet another serendipity is that the divergence of the 
bare potential is only logarithmic, cf. Eq. (8). The en- 
ables us to compute the non-perturbative long-range re- 
sponse with a logarithmic accuracy. Indeed, in view of 
Eq. (10), $k and Vk are linearly related at all k, ex- 
cept perhaps k <~ l/R. But this "difficult" range of in- 
termediate k makes an C(£" 1 ) contribution to the to- 
tal l>(r). To this order in £ we can treat <l(r) as a 
Gaussian random potential. Its probability distribution 
function (PDF) is uniquely determined by the correla- 
tor K{r) = ($(0)$(r)), to calculate which we can use 
Eq. (10). Since K(r) ~ K(r) for rA > 1, this immedi- 
ately leads to Eq. (5). 

Next, Gaussian statistics implies ergodicity. Therefore, 
the density correlation function can be written as 

S{r) - J dfdf P M (r, /, /' )«(/) n(f ) , (15) 

where = (S(f - <&(0)) S(f - #(r))) is the two-point 
PDF of <f> and function n(f) is given by Eq. (13). Using 
the standard expression [27] for in terms of K and a 
bit of algebra [23], I get Eq. (3). At r » R, it simplifies 
to S(r) ~ (x) 2 K(r), where ( X ) = J dfP*(f) X (f), which 
implies that xo = (x)- Substituting here Eq. (14) and 
computing the average over the Gaussian field $ using 
Eq. (5), I get Xo = 8if(^)/(7r 3 ftV). Finally, since K(l) 
is given by Eq. (4), we recover the desired Eq. (2). 

It is also possible to compute the PDF P n of n. The 
Gaussian statistics of $ combined with Eq. (13) entails 



P n {n) = 



2 V / 7r£T 



exp 



\n\f 
£ 



|n|^ 2 >l. (16) 



Note, however, that this equation is invalid at small den- 
sities, |n| < l/£ 2 , where the RTFA fails. At such n the 
divergence of P n (n) basically saturates [23]. These low- 
density regions arc usually found near the PNJ. 

Geometrically, the PNJ are the isolincs $(r) — of a 
surface with the "height" profile $(r). The zero height 
is the percolation threshold, and so all but one of the 
PNJ are closed loops. These loops are characterized by a 
certain fractal dimension D^. For loops of diameter d in 
the range A -1 <C d <SC R, in which <l(r) surface is loga- 
rithmically rough [Eq. (5)], we have the exact result [28] 
Dh = 3/2, which means that closed- loop PNJ typically 
have the perimeter length of 

P ~d 3/2 r 1/2 , A- 1 <d<p. (17) 
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As d increases beyond R, the correlator K(d) rapidly 
decays [Eq. (5)], and so Dh crosses over [28] to the usual 
uncorrelated percolation exponent [27] of 7/4. 

Let us now discuss electron transport. Away from the 
NP where electron density n is large and homogeneous, 
transport can be studied by means of the usual kinetic 
equation. The final results for conductivity cr(n) and 
transport mean-free path l(n) are [6, 8, 12, 13] 

a(n)=-^£ 2 \n\, l(n) = -^ £ 2 \n\ 1 / 2 . (18) 

7T h 7T ' 

Our goal is to compute the conductivity ctnp at the NP 
where n = n(r) is inhomogencous. 

The first step is to show that we can define the con- 
ductivity a(r) and the mean- free path Z(r) locally. Un- 
der the assumed condition 1/R <C A < l/l for typi- 
cal n(r) <~ C/£ 2 we have Vl(h) ~ Al(h) < 1. Hence, 
function l(h(r)) is slowly varying and Eq. (18) for the 
uniform density can be used: I = 1(h). In turn, 
the local conductivity is given by the Einstein relation 
<r(r) = e 2 x(r)vl(r)/2. To be careful one should check 
this zeroth order result by calculating higher order cor- 
rections due to spatial fluctuations of the collision term 
in the kinetic equation. This can be done by treating 
<5Z — 1 (r) = £ _1 (n(r)) — l _1 (h) as a perturbation. Skip- 
ping the details [23], I just announce the conclusion: if 
\h\ 3> l/l 2 , then the fluctuations of <5Z _1 (r) are self- 
averaging on the scale of 1(h). In other words, correc- 
tions to 1(h) are parametrically small. For simplicity, I 
will ignore them together with the similar corrections to 
x(r) [the second term in Eq. (14)] to obtain 

&(v) — o~(h(v)) , \h(v)\-»l/l 2 . (19) 

Problems arise at \h\ < l/£ 2 where the corrections to 
Eq. (19) exceed a(h). In the classical regime a ^> e 2 /h, 
we could have handled this by adopting a model form 

a(r) = aoC(r) , \h\ < l/£ 2 , (20) 

where C(r) > is some random function with the corre- 
lation length A -1 ~ I and the typical value of the order 
of unity [29]. However, in graphene cto = a(\j£ 2 ) — 
(8/w)e 2 /h is so low that the very concept of local con- 
ductivity is potentially jeopardized by quantum interfer- 
ence and localization effects. Fortunately, Eq. (20) is 
saved by the special geometry of the regions where it 
is intended to be used. These regions typically form 
ribbons of width x <~ £ that follow the PNJ loops. 
These low-conductance "ribbons" are connected to high- 
conductance "reservoirs" on both sides of the PNJ. In 
such a geometry the importance of the quantum effects 
is controlled [32] by the total conductance G p across the 
perimeter length p of the PNJ. For p ^> £, it is given by 

G p ~<7of»^, (21) 



in which case the scmiclassical model (20) is justified. 
Again we have succeeded in reducing the complicated 
quantum problem to a simpler classical one: finding 
the macroscopic conductivity cnp of an inhomogencous 
medium with local conductivity given by Eqs. (19) and 
(20) as a function of the local density h(r). 

A rigorous upper bound on ctnp is the spatially aver- 
aged conductivity. Due to ergodicity of Gaussian fields, 
the averaging can be done over h instead. Using Eqs. (16) 
and (18), I obtain: 

cr N p < (5-) = J dhP n (h)a(h) ~ (4£/n)e 2 /h. (22) 

(Incidentally, this enables me to conclude that Ref. 13 
overestimates ctnp in the limit a « 1.) 

Now I present the argument crucial to my theory of 
transport. It shows that ctnp is not sensitive to the de- 
tails of the model form (20). Instead, it is determined by 
the typical local conductivity, and so is not far below the 
upper bound (22). This statement should be contrasted 
with other theoretical views on the subject. Since the 
low-conductivity regions (20) reside at the percolation 
contour, it has been suggested [1, 13, 33] that the trans- 
port at the NP may be governed by percolation [27]. In 
that picture the current paths are severely constrained in 
order to avoid crossing the PNJ as much as possible and 
ctnp depends on their average transparency [33]. 

I show that in the model under study the closed loop 
PNJ contours are not resistive enough for the percolation 
effects to develop. Indeed, the percolation approach [33] 
would apply only if there existed a wide range of loop 
diameters d ~> I such that the conductance G p of the 
loop perimeter were much lower than the conductance 
Gd of their interior. Suppose first that d <C R for such 
loops. Let us show that it leads to a contradiction. If 
the current avoids crossing the perimeter, it has to flow 
in narrow channels of some width w <C d due to the 
fractal geometry of the loop. Hence, Gd ~ a-^w/d <C 
<7i n t where (Tint is the typical local conductivity inside the 
loop. To estimate the latter, I note that from the general 
properties of Gaussian random fields [27] the potential <l 
at a distance r from the PNJ (l> = 0) has the variance 

<* j >™^>-^K(t) 2 >"G). m 

so that the corresponding density is h ~ ln(r/£)/£ 2 . Us- 
ing r ~ d and Eq. (18), I get cri nt ~ a ln(d/£). In com- 
parison, G p - o-v(d/£f/ 2 , cf. Eqs. (17) and (21). Hence, 
G p 3> Gd, which contradicts the assumption made (the 
same is true for d 3> R). I conclude that in the present 
model the percolation-type transport is not realized. 

In the absence of strong ramification of current paths 
by the PNJ, I expect the numerical coefficient c in Eq. (7) 
to be of the order of unity. To estimate it more accu- 
rately I use the effective medium theory (EMT). Indeed, 
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the EMT is usually adequate for systems where percola- 
tion effects are unimportant. Within the EMT, ctnp is 
determined from a certian nonlinear equation, the two 
most popular versions of which were originally proposed 
in Ref. 34 and Ref. 35: 



ctnp 



a + (D — 1)ct np 



/ 



dz I 
— In ( cxp 



V ctnp 



1 

V 



V 



(Bruggcman) , (24) 



(Hori) . 



(25) 



Here V = 2 is the space dimension. Equations (24) 
and (25) can be viewed [35] as an approximate resum- 
mation of the infinite diagrammatic series for the macro- 
scopic conductivity using, respectively, the self-consistent 
single-site approximation and the cumulant expansion. 

The averages in Eqs. (24) and (25) are dominated by 
typical ct, so it suffices to use P n (n) from Eqs. (16). 
A straightforward numerical solution of these equa- 
tions then gives: CTNp(Bruggeman) = 0.48e 2 //i and 
ctnp (Hori) = 0.55e 2 //i. Taking the difference of the two 
as a measure of their accuracy, I arrive at Eq. (7) . 

Finally, I briefly comment on experimental implica- 
tions of the presented theory. For simplicity, I have as- 
sumed the dielectric constant of the medium to be the 
same on both sides of graphenc. It is more realistic to 
have graphene at the interface of a half-space with dil- 
electric constant k\ <~ 1 and a film of dielectric constant 
k 2 ^> 1 and thickness D, whose other side is covered 
by a metallic gate. In this geometry k is replaced by 
(ki + K2)/2 while C becomes min( ln(R/£) , ln(D/£)). 
Small a can be achieved experimentally using ice [4], 
ethanol, and other dielectrics [36]. However, it is dif- 
ficult to make C larger than C ~ 5, for which Eq. (7) 
predicts ctnp ~ 2.5e 2 /h. This is two to three times lower 
than the measured [1, 2, 3, 4] ctnp- This may indicate 
that in experiment the charged impurities arc cither not 
exactly coplanar with graphene, or are correlated, or are 
not the only source of disorder. The last possibility is 
further corroborated by a very modest increase in mobil- 
ity [36] away from the NP upon a large in situ increase 
of K. 
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► Equation (12). This equation can be written as 



3e 



3 4 m 



T 



(12) 



The first term is the energy density of the state with uniform density n. All interaction effects are assumed to be 
included in the renormalized Fermi velocity v. In principle, this makes v a function of n. However, this function is 
very slow [1], and so we treat it as a constant. The correction term 5e arises due to interaction of electrons with 
large-Zc harmonics N k of the impurity charge density. The corresponding energy density can be expressed in terms of 
the linear-response dielectric function e(k) of graphene: 



Ss = 



1 

2^2 



£ 

fe>A 



e(k) 



C4 



e(oo) 



(N k N_ k ) = 



hv 

8q^ 



dk 



1 

e(fc) 



1 



e(oo) 



(SI) 



The subtraction in the integrand ensures convergence. (It amounts to a constant shift of no significance for the 
problem studied.) Suppose that A ^> Aakp, then the dominant contribution to the integral comes from momenta 
A < k < kp, at which e(k) ~ 1 + {Aakp/k). Using this approximation and replacing the upper integration limit by 
fcp, we obtain Eq. (12). A more careful calculation that utilizes the full expression [2] for e(k) confirms this result. 
► Equations (3), (15), and (23). In the main text I allude to a "standard expression" for the two-point probability 
distribution function (PDF) of i>(r). This expression is given by [3, 4]: 



exp 



K{r) K(0) 2 2 

Q(r) JJ 2Q(r) {J J ' 



, Q[r)=K\r)-K\Q) 



(S2) 



It is easy to see that it respects the identity (<I>(0)<I>(r)) = K (r) and gives the correct PDF of a local <!>(r): 



P*(f) = j <V'P**(r,f,f') = 



exp 



V 



(S3) 



On the other hand, if the variance of <l is computed using the constrained PDF P$$(r, /, /' = 0), Eq. (23) is obtained. 

Next, substituting Eq. (S2) into Eq. (15) of the main text and directly evaluating the double integral, I get Eq. (3). 
Finally, to derive the relation between S(r) and K(r) at large r, I use the expansion 



P**(rJJ')-P*(f)P*(f')-K(r) 



dP*(f) dP*(/') 



df 



df 



r > R. 



(S4) 



This formula is rather general and is not limited to Gaussian random fields. Indeed, it can be derived without relying 
on Eq. (S2) but using the primary definition of the two-point PDF: = (5 (/ — $(0)) <5(/' — $(r)) . The idea is that 
the statistical dependence of random variables $(0) and l>(r) is due to Fourier harmonics <&k with kr < 1. At large r 
such <J>k have a small phase space and can be considered a perturbation to a local <!. Therefore, one can expand the 
arguments of the (5-functions to the second order in such <!>k and average the result to get Eq. (S4) . 
Substituting Eq. (S4) into Eq. (15) and integrating by parts twice, I get 



S{r)~K{r) / dftf P*(f)P*(f 



dn(f) dn(f') 
df df 



(x) 2 K(r), 



(S5) 



2 



from which the key result xa = (x) follows. Using Eqs. (S5) and (14) on can show that xo does not depend on the 
arbitrary cutoff A and is given by Eq. (2) of the main text. 
► Equation (16). This equation reads 

Pn(n) = —!= exp f-^n , \n\f » 1 . (16) 



n 



It directly follows from Eq. (S3) and the renormalized Thomas-Fermi approximation (RTFA), Eq. (13), if the second 
(correction) term therein is neglected. Here I wish to discuss how Eq. (16) may be modified at |n| < l/£ 2 where the 
RTFA fails. According to Eq. (16), such h are much smaller than the typical n ~ C/£ 2 . They are usually found near 
the p-n junctions (PNJ), i.e., near the contours 3>(r) = 0. This is because by symmetry Vn and n are independent 
random variables, and so |Vn| is generally not small when n is. Therefore, for small n the distance x — n/|Vn| to the 
nearest PNJ is also small. Equation (23) implies that in this case l>(x) = (hv/l) \Jir\n (xA) £, where £ is a Gaussian 
random variable with zero mean and unit variance. The RTFA, which is a local theory, fails for x smaller than the 
Fermi wavelength l/kp(x) ~ hvj |$(x)|. This failure occurs at x < X£ ~ I / assuming of course that x^ <C R, 
i.e., that £, which is specific for each PNJ, is not too small. Based on our previous study of graphene PNJ [5, 6], we 
expect that at x < x^, where \n\ < \n(x^)\ <~ i 2 /£ 2 , the density profile is linear: |n(x)| ~ x/x|. The local inverse 
compressibility in the same domain is given by <~ hvx^. The corresponding P n {n) can therefore be written as 

P„(n) * P $ (0)((x~ 1 )) ~ e 2 £~ 1/2 « » ~ ^= In ^ , ^ « N « j 2 . (S6) 

Here ((. . .)) stands for the averaging over £ under the indicated above condition £ 2 > \n\£ 2 . Equations (16) and (S6) 
match at their common border of validity \n\£ 2 ~ 1. At |n| <C 1/-R 2 , the divergence of P n (n) presumably saturates. 
► Equation (19). This equation states that the local conductivity is determined by the mean-free path I computed 
using a coarse-grained density h, which is, in turn, obtained by averaging local density n(r) over a neighborhood 
of a certain size A" 1 . Here I show that the proper A" 1 is of the order of the typical mean- free path £y/~C. To this 
end consider a region of the sample away from the PNJ where n(r) is large enough and the semiclassical kinetic 
equation vVf((f), r) = St[/] + cos cf) applies. Here / is the deviation of the quasiparticle distribution function from the 
equilibrium, <j) is the angle between the directions of quasiparticle velocity v and the electric field, and St[/] is the 
collision term. The importance of density fluctuations can be estimated using a simplifying assumption of isotropic 
scattering, St[/] = — ^ 1 (r)/, and then treating <5Z _1 (r) = Z _1 (n(r)) — l^ 1 (n) as a perturbation. Applying the method 
of characteristics, I get the leading-order result 



1 = JdrG(r)(cxp - J dl^iX^rdX 



_3/2 1 



Expanded to the the order O^r 1 ), it reads I = l(n) + l 2 {n) J dr g(r)(5l- 1 (0)5r 1 (r)}, which can also be derived 
diagrammatically. It is easy to see then that the difference between I and l(n) is small when the relative change in 
n(r) over the distance of l(n) is much smaller than unity. In turn, this is satisfied when the same can be said about 
the potential <&(r). On the other hand, Eqs. (S2) and (13) entail 



$(r) - $(0) 



$(r) 



'Q(r) ^ I ln(r /I) (g8) 



|*(r)| V K(0) V 



For r = l(n) [Eq. (S7)] the right-hand side of this expression is indeed much smaller than unity as long as \n\£ 2 >• 1, 
which is precisely the domain of validity of Eq. (19) indicated in the main text. 
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